clear
set more off

local critbase "& age >= 30 & age <= 45"
local crittrunc "& age >= 30 & age <= 40"
local critextend "& age >= 30 & age <= 45"
local critcomb "& age >= 30 & age <= 45"

local yearbase $yearlist
local yeartrunc $yearlist
local yearextend $yearlist
local yearcomb $yearlistcomb

local varlist1 $varlist1
local varlist2 $varlistcomb


local yearlist $yearlist
local varlist1 $varlist1
	
*************************************************
*Individual Level Regressions
*************************************************

*Get data

use $tempdir/panel_trends_30yr_exage.dta, clear

keep if race==1
keep serial pernum statefip datanum first year1 year2 statefip bpl spousebpl urban age spouseage agegap sex race spouserace married ever_married logocc* qtl* husboccscore spousefirst

gen logoccscore_father_sp_base = logoccscore_father_spouse if spouseage >= 30 & spouseage <= 45
gen logoccscore_father_sp_trunc = logoccscore_father_spouse  if spouseage >= 35 & spouseage <= 45
gen logoccscore_father_sp_extend = logoccscore_father_spouse if spouseage >= 30 & spouseage <= 55

gen logoccscore_father_sp_comb = logoccscore_father_spouse if (spouseage >= 30 & spouseage < 40 ) | (spouseage >= 40 & spouseage <= 45 & logoccscore_father_spold == .)
replace logoccscore_father_sp_comb = (logoccscore_father_spouse + logoccscore_father_spold) / 2 if spouseage >= 40 & spouseage <= 45 & logoccscore_father_sp_comb == .
replace logoccscore_father_sp_comb = logoccscore_father_spold if spouseage > 40 & spouseage <= 55 & logoccscore_father_sp_comb == .
replace logoccscore_father_sp_comb = . if inlist(year2,1880,1930)

order year2 sex spousefirst spouseage spouseage logoccscore_father_spouse logoccscore_father_spold logoccscore_father_sp_extend logoccscore_father_sp_comb
sort year2 sex spousefirst spouseage

*Variables for the relevant marriagemarket and childhood years
gen maryear = year2 - 10
gen chyear = year2 - 30

*Here, assign the state the match the panel controls as statefip, state of residence
*Change to bpl to match on birth place instead
gen fips = statefip

* Regressions of ever married status
egen firstcode = group(sex first year2)

gen int bpl2 = bpl
replace bpl2 = 9999 if bpl>100

gen byte foreign_born = bpl2>100 if bpl2~=.

sort year2 foreign_born

gen urb = urban==2 if urban!=.
label var urb "Urban dummy"


forvalues i = 1/4 {
	gen q`i' = qtl_score_wgt_full==`i'
	replace q`i'=. if qtl_score_wgt_full==.
}

foreach y in `yearbase' {
		gen y`y' = year2==`y'
}

forvalues i=1/4 {
	foreach y in `yearbase' {
		gen q`i'_`y' = q`i'*y`y'
	}
}

gen logoccscore_husband_temp = log(husboccscore)

local dep_ever_married = 1
local dep_logoccscore_dadinlaw_temp = 2

local dep_base = 1
local dep_trunc = 2
local dep_extend = 3
local dep_comb = 4

gen logoccscore_dadinlaw_temp = . 

cap postclose results
postfile results agespec depvar spec qtile year b se using "$tempdir/reg_robust_exage.dta", replace

foreach spec in base trunc extend{

replace logoccscore_dadinlaw_temp  = logoccscore_father_sp_`spec'

foreach var of varlist ever_married logoccscore_dadinlaw_temp {

* Ever married regressions -- no bpl controls
reg `var' `varlist1' i.age foreign_born if sex==2 `crit`spec'', cluster(firstcode)
gen s1_`var'`dep_`spec'' = e(sample)

forvalues i=2/4 {
*	local b`i'_1880 = _b[q`i']
	foreach y in `year`spec'' {
		local b`i'_`y' = _b[q`i'_`y']
		local se`i'_`y' = _se[q`i'_`y']
	}
}

forvalues i=2/4 {
	foreach y in `year`spec''  {
		post results (`dep_`spec'') (`dep_`var'') (1) (`i') (`y') (`b`i'_`y'') (`se`i'_`y'')
	}
}


* Ever married regressions -- full bpl controls
reg `var' `varlist1' i.age i.bpl2##i.year2 if sex==2 `crit`spec'', cluster(firstcode)
gen s2_`var'`dep_`spec'' = e(sample)

forvalues i=2/4 {
*	local b`i'_1880 = _b[q`i']
	foreach y in `year`spec''  {
		local b`i'_`y' = _b[q`i'_`y']
		local se`i'_`y' = _se[q`i'_`y']
	}
}

forvalues i=2/4 {
	foreach y in `year`spec''  {
		post results (`dep_`spec'') (`dep_`var'') (2) (`i') (`y') (`b`i'_`y'') (`se`i'_`y'')
	}
}

}
}

foreach spec in comb {

replace logoccscore_dadinlaw_temp  = logoccscore_father_sp_`spec'

foreach var of varlist ever_married logoccscore_dadinlaw_temp {

* Ever married regressions -- no bpl controls
reg `var' `varlist2' i.age foreign_born if sex==2 `crit`spec'', cluster(firstcode)
gen s1_`var'`dep_`spec'' = e(sample)

matrix vcv = e(V)

forvalues i=2/4 {
*	local b`i'_1900 = _b[q`i']
	foreach y in `year`spec'' {
		local b`i'_`y' = _b[q`i'_`y']
		local se`i'_`y' = _se[q`i'_`y']
	}
}

forvalues i=2/4 {
	foreach y in `year`spec'' {
		post results (`dep_`spec'') (`dep_`var'') (1) (`i') (`y') (`b`i'_`y'') (`se`i'_`y'')
	}
}


* Ever married regressions -- full bpl controls
reg `var' `varlist2' i.age i.bpl2##i.year2 if sex==2 `crit`spec'', cluster(firstcode)
gen s2_`var'`dep_`spec'' = e(sample)

matrix vcv = e(V)

forvalues i=2/4 {
*	local b`i'_1900 = _b[q`i']
	foreach y in `year`spec''  {
		local b`i'_`y' = _b[q`i'_`y']
		local se`i'_`y' = _se[q`i'_`y']
	}
}

forvalues i=2/4 {	
	foreach y in `year`spec''  {
		post results (`dep_`spec'') (`dep_`var'') (2) (`i') (`y') (`b`i'_`y'') (`se`i'_`y'')
	}
}



}

}

postclose results

save "${tempdir}reg_data_exage.dta", replace 

clear 
save "${tempdir}reg_nobs_exage.dta", replace empty

local dep_ever_married = 1
local dep_logoccscore_dadinlaw_temp = 2
local dep_logoccscore_husband_temp = 3
local dep_agegap = 4

foreach x in 1 2 3 4{
foreach var in ever_married logoccscore_dadinlaw_temp{
foreach spec in 1 2 {

use "${tempdir}reg_data_exage.dta", clear

keep if s`spec'_`var'`x' == 1

sort year2 first
by year2 first: gen nnames_female = _n == 1

if "`var'" == "logoccscore_dadinlaw_temp" {
sort year2 spousefirst
by year2 spousefirst: gen nnames_male = _n == 1
	}
	
collapse (sum) nobs=s`spec'_`var'`x' nnames_*, by(year2)	

gen agespec = `x' 
gen depvar = `dep_`var''
gen spec = `spec'

rename year2 year 

append using "${tempdir}reg_nobs_exage.dta"
save "${tempdir}reg_nobs_exage.dta", replace 

}
}
}

use "$tempdir/reg_robust_exage.dta", clear
merge m:1 agespec depvar spec year using "${tempdir}reg_nobs_exage.dta"

export excel using "$outputdir/Regression Robustness.xlsx", sheet("exage", replace) firstrow(var)

erase "${tempdir}reg_data_exage.dta"
erase "${tempdir}reg_nobs_exage.dta"




